!=========Programa que calcula a seção de Poincaré do pêndulo===========

PROGRAM exerE

!=======================================================================

IMPLICIT NONE

!=======================Declaração de variáveis=========================

      INTEGER :: I, J, steps, flag
      REAL(8) :: tethai, tethai1, omegai, omegai1, Dt, currt
      REAL(8) :: l, totalt, Omega, gam, F0, tetha0
      REAL(8), PARAMETER :: g = 9.80665d0, PI = 3.14159265358979d0

!=============Inicialização, leitura e abertura de arquivos=============

      WRITE(*,*) "Entre com o comprimento da haste:"
      READ(*,*) l
      WRITE(*,*) "Entre com o intervalo entre iterações delta_t:"
      READ(*,*) Dt
      WRITE(*,*) "Entre com o ângulo inicial do pêndulo (em graus):"
      READ(*,*) tetha0
      WRITE(*,*) "Entre com o tempo total de simulação:"
      READ(*,*) totalt
      WRITE(*,*) "Entre com o parametro resistivo gamma:"
      READ(*,*) gam
      WRITE(*,*) "Entre com a amplitude da força externa F_0:"
      READ(*,*) F0
      WRITE(*,*) "Entre com a frequência da força externa:"
      READ(*,*) Omega

      !Abre arquivo para omega(tetha)
      OPEN(UNIT=10, FILE='exerE_out.dat', STATUS='unknown')

      !Converte o angulo de graus para radianos
      tetha0 = tetha0*PI/180.0d0

      tethai1 = 0.0d0
      omegai = 0.0d0
      omegai1 = 0.0d0
      currt = 0.0d0
      flag = 0
      J = 0

!=============================Faz as iterações==========================

      steps = totalt/Dt
      DO
            currt = 0.0d0
            tethai = tetha0
            omegai = 0.0d0
            flag = 1
!            write(*,*) "teste"
            DO I=1,steps
                  !Escreve no arquivo
                  IF (ABS(currt - 2*J*PI / Omega) < (Dt/2.0d0)) THEN
                        WRITE(10,*) tethai, omegai
                        !Muda a flag para iterar o proximo J
                        flag = 0
                  END IF
                  !Faz a iteração de omega
                  omegai1 = omegai - (g/l)*DSIN(tethai)*Dt - gam*omegai*Dt +&
                  &F0*DSIN(Omega*currt)*Dt
                  !Faz a iteração de tetha
                  tethai1 = tethai + omegai1*Dt
                  !Muda os ângulos para deixar entre -pi e pi
                  IF (tethai1 .GE. PI) THEN
                        tethai1 = tethai1 - 2.0d0*PI
                  ELSE IF (tethai1 .LE. -1.0d0*PI) THEN
                        tethai1 = tethai1 + 2.0d0*PI
                  END IF
                  !Calcula em qual tempo ocorreu essa iteração
                  currt = currt + Dt
                  !Prepara variavaies para a proxima iteração
                  omegai = omegai1
                  tethai = tethai1
            END DO
            J = J + 1
            IF (flag .EQ. 1) EXIT
      END DO

!=======================================================================

END PROGRAM exerE

!=======================================================================
